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1 . INTRODUCTION 


Del ami nati on is a mode of failure that is unique to composite 
laminates. A delamination crack can significantly reduce the 
compressive strength of a laminate, and can also induce matrix cracking, 
which will further degrade the structural integrity of the laminate. 
Delamination can be produced by both static and dynamic loads. Great 
attention has been given to the problem of free-edge delamination in 
laminates subjected to in-plane static and fatigue loadings [1-4], and 
many people have attempted to measure the fracture toughness with 
respect to delamination cracks [5-9]. Few attempts have been made, 
however, to measure the dynamic delamination fracture toughness. This 
is one of the goals of this research. 

It is known that material properties often exhibit a strain-rate 
dependence [10]. In fact, the fracture toughness measured in a static 
loading test will, in general, be different from that observed under 
dynamic loading [11, 12]. Use of a statically measured fracture 
toughness in a dynamic crack propagation analysis can indeed seriously 
overestimate the crack arrest capability of a structure [13]. In 
addition, the effective initiation toughness measured under impact 
loading conditions [14] was found to be roughly twice that measured in 
static loading tests of high strength steel. The differences between 
the static and dynamic toughness are due primarily to the fact that the 
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static calculations cannot account for the return of kinetic energy to 
the crack tip, and the assumption that the fracture energy is indepen- 
* dent of the crack propagation velocity [15]. Early work in crack 
arrest and dynamic fracture mechanics in general was motivated by 
brittle fracture commonly observed in the hulls of ships [16, 17]. To 
determine if the arrest toughness, K_, was a material characteristic, 

a 

Hoagland [16] performed experiments on wedge-loaded double-cantilever 
beam specimens of four different steels. The amount of stored strain 
energy at the onset of crack propagation was systematically varied 
from one specimen to another by changing the radius of the starter 
notch. As the initial notch tip bluntness was increased, the stress 
intensity at initiation increased, and a resulting decrease in stress 
intensity at arrest was observed. From these results it was concluded 
that the "arrest toughness" is not a material constant, but is 
dependent on the amount of initial strain energy stored in the 
structure. 

An energy balance criterion can account for the return of kinetic 
energy to the crack tip. Hahn [17] concluded that fast fracture and 
arrest in a variety of steels are governed by an energy balance 
criterion, in which the excess energy (G-R) is available to drive the 
crack tip. When G<R, the kinetic energy contributes significantly to 
the crack driving force until arrest. In [18] Hahn studied fast 
fracture in a wedge-loaded double cantilever beam specimen using a 
beam-on-elastic-foundation model. He showed that there is a unique 
relationship between the steady-state crack propagation velocity and 
the crack length at arrest for a given material and specimen geometry. 
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Indeed, the stress intensity at arrest was found to vary with crack 
propagation history, and therefore it cannot be considered a material 
* property. 


A review of dynamic fracture toughness measurements in polymers 
is given by Kobayashi [19]. He points out that evidence from 
experiments on large polymeric specimens [20-22] seems to suggest a 
unique relation between the dynamic fracture toughness, K T , for a 
propagating crack, and the crack propagation velocity, in polymers. 
Hodulak [23] tested dynamic fracture specimens of three different 
geometries to assess the uniqueness of the fracture toughness versus 
crack velocity relationship. He used a finite element analysis with 


an experimentally determined K. vs. a relationship, and calculated the 

D 

crack propagation history. His results showed that K T vs. a relation 

D 

is geometry dependent, and he concluded that dynamic crack propagation 
is not controlled solely by the instantaneous state surrounding a 
running crack, but is also significantly affected by the motion of the 


structure remote from the crack. 


In contrast to in-plane static loads, under which delamination 
often initiates from free edges, impact loading always results in 
interior delamination near the impact zone. Thus, the delamination 
mechanism cannot be explained by using the free edge singular stress 
concept. The strain energy release rate was defined by Erdogan [24]. 
Several methods of calculating the energy release rate from a numerical 
analysis have been proposed. Rybicki [25] devised an efficient technique 
to calculate static Mode I and II energy release rates from a single 
finite element analysis. He used the near-tip nodal forces and 
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displacements to calculate the work required to close the crack by an 
amount Aa, assumed equal to one near-tip element length. It has been 
shown [26] that this crack closure technique is also valid in the case 
of dynamic loading, if the finite element mesh near the crack tip is 
small enough. The first objective of this work is to model the impact 
of a composite laminate and to determine from this model a critical 
value of strain energy release rate, G c , required to cause instability 
of an existing delamination crack. This parameter, if it exists, may 
be a characteristic of the material. To model the impact response of 
the laminate adequately, an accurate representation of the impact force 
history must be determined. This is the second major objective of 
this work. 

Insofar as the primary concern here is to model the dynamic 
behavior of the laminate from impact to the initiation of propagation 
of the delamination crack, no attempt is made to model the propagation 
of the crack. Therefore, the analysis of the cracked laminates 
presented here is valid only until the onset of crack propagation in 
the laminate. 

The presentation of this work is organized as follows. Chapter 2 
explains the equipment and procedures used to perform the Ballistic 
Impact experiments. In Chapter 3, the photographic data and the 
corresponding measurements taken from it are presented. A technique 
used to characterize the stiffness loss in the delaminated impact 
specimens is illustrated. Chapter 4 describes how the impact force 
history was modeled, and how this model was incorporated into the 
impact analysis of the composite laminate. Chapter 5 is a presentation 
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and discussion of the results of the impact analysis of the laminates. 
Chapter 6 is the summary. 
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2. EXPERIMENTAL APPARATUS AND PROCEDURE FOR IMPACT TESTING 


2.1 Impact Specimen Preparation 

Impact specimens were cut from 20-ply [90/0]g s T-300/934 graphite/ 
epoxy laminates of dimension 0.105x12x18 in. These laminates were 
fabricated at the NASA Lewis - Research Center in Cleveland, Ohio. A 
delamination crack was embedded in each laminate by placing a 
.001x1.0x18 in. strip of Teflon between two plies during the lay-up 
process. This prevented the two adjacent plies from bonding together 
in this area. A beam- like geometry was chosen for the impact 
specimens. Nominal dimensions are shown in Fig. 2.1. Thus, the 
initial delamination is a 1.0 inch long, through-the-width crack. The 
location of the embedded crack in both the longitudinal and thickness 
directions was varied between laminates. This was done to study the 
effect of crack location on delamination characteristics. 

2.2 Ballistic Impact Testing of Composite Specimens 

Silicon rubber balls i inch in diameter were used as impactors. 
These relatively soft impactors do not cause significant surface 
damage near the impact site, thus allowing crack extension to be the 
primary mode of impact damage. Nitrogen gas was used to fire the 
impactor through the cannon, shown in Fig. 2.2. A chamber pressure of 
20 psi could propel the 1 gram rubber ball at approximately 6000 inches 
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per second. The impact velocity was determined by two pairs of photo- 
electric diodes placed on both sides of the path of the impactor near 
* the muzzle of the barrel. A high speed 16 mm FASTAX framing camera 
was used to record the crack propagation. It was mounted to give an 
edge-on view of the impact specimen, which was enclosed in a plexiglas 
box and clamped at one end in a cantilever fashion. The peak framing 
rate of the camera is 8000 frames per second. This rate is effectively 
doubled by an internal rotating prism which made two exposures per 
frame, thus taking 16000 pictures per second. Because of the high 
exposure rate of the film, very bright light was needed to adequately 
illuminate the impact specimen. This was provided by three 100-watt 
floodlights. The firing sequence was initiated from a control panel 
with timers set to trigger the camera and photo lights just before 
impact. 

The details of this experimental set up are shown in Figs. 2.3 
and 2.4 

2.3 Measurement of Delamination 

To facilitate the measurement of the delamination length from the 
high speed film, the edge of the impact specimen nearest the camera 
was painted yellow. This made the advancing crack clearly visible 
against the light background. In addition, dark stripes were painted 
at quarter-inch intervals along the edge of the specimens to serve as 
reference points in measuring the crack length. The Ultrasonic C-Scans 
of several impacted specimens shown in Fig. 2.5 indicate that the 
variation in delamination crack length through the width of the 
specimen is small. The measurements of crack length taken along the 
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outer edge of the specimen can therefore be assumed to represent the 
total delaminated area at a given time. 




Figure 2.1 Nominal Impact Specimen Dimensions 
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Figure 2.2 Ballistic Impact Experimental Set-Up 
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Figure 2.3 Camera Arrangement 
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Figure 2.5 Ultrasonic C-Scans of Impacted Specimens 
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3. IMPACT RESPONSE OF CRACKED LAMINATES - 
EXPERIMENTAL RESULTS 


3.1 Threshold Impact Velocity 

The dependence of delamination damage on impact velocity is of 
primary interest. In this study, attention will be focused on the 
threshold velocity at which the embedded crack becomes unstable. 

Figure 3.1 shows the geometry of eight different specimen config- 
urations tested. The location of the impact point varied slightly 
between specimens due to small variations in alignment of the gun 
barrel. The extent of this variation, along with the relation between 
impact velocity and total delaminated area for each of the specimen 
configurations tested is shown in Tables 3. 1-3.7 and in Figs. 3. 2-3. 7. 
Each specimen contains an initial (embedded) delamination. In most 
cases, the existence of an unambiguous threshold velocity is quite 
evident. Threshold velocities for specimen configurations A-H were 
determined from Figs. 3. 2-3. 7. Insufficient data was available for 
specimen configurations F and G so corresponding threshold velocity 
plots are not shown for these configurations. Among the three 
thickness locations tested, threshold velocity is greatest for the 
mid-plane crack (Table 3.2, Fig. 3.3) and lowest for the lower 
off-midplane crack (Table 3.3, Fig. 3.4). The results show that when 
impacted near the crack tip, the delamination crack becomes unstable 


14 



at lower velocities. Tables 3.4 and 3.5 show that this phenomenon is 
more pronounced for cracks located near the top (impact) surface. This 
'behavior may be a result of the unsymmetrical distribution of shear 
stress over the cross section of the laminate which occurs near the 
point of impact. Joshi [27] showed that the maximum shear stress 
occurs near the impacted side of the laminate cross section during and 
shortly after the contact interval. This would suggest that the onset 
of crack propagation in specimens of the geometry shown in configura- 
tion D of Fig. 3.1 is dominated by a shearing (Mode II) rather than an 
opening (Mode I) action. 

3.2 Dynamic Crack Propagation 

3.2.1 Midplane Delamination 

A typical impact sequence is shown in Fig. 3.7. Characteristics 
such as duration of contact period and beam displacement response can 
be estimated from the figure. It should be noted that all measurements 
were taken from larger images projected on a movie screen. The figures 
shown here are primarily for illustration. In this case, the embedded 
crack lies along the specimen midplane and directly under the impact 
site, as shown in the figure. The resulting crack propagation is shown 
in Fig. 3.8. The crack arrest (438 < t < 688ys) is apparently due to 
the nature of strain response near the propagating crack tip. A 
decrease in local curvature of the beam is accompanied by a decrease 
in available crack driving force. This correspondence is shown in 
frames 9-11 of Fig. 3.7. Frames 12-14 (688 < t < 813us) show the 
subsequent increase in curvature, and the corresponding resumption of 
crack propagation. 
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Apparently, the geometry of the impact specimen can significantly 
affect the crack propagation. Strain (curvature) will be affected by 
the arrival, of flexural wave reflections from the boundaries, so the 
position of the crack relative to the boundaries will affect crack 
propagation. The time delay between impact and initial crack propaga- 
tion observed in Figs. 3.7 and 3.8 is a result of the impact occurring 
directly on the embedded crack. The distributed compression on the 
crack faces caused by the deforming impactor (63 < t < 375us, Fig. 3.7) 
prevents any crack propagation from occuring during the contact 
interval. Now, if the embedded crack is moved sufficiently away from 
the impact site, as depicted in Figs. 3.9 and 3.10, the interference 
of the impactor with crack propagation should be minimized. Compare 
Figs. 3.7 and 3.8 with Fig. 3.9. All three figures show similar crack 
arrest characteristics as the reflections arrive. However, Figs. 3.9 
and 3.10 show a significant difference in time between impact and 
onset of crack propagation. 

Impact specimens with different initial crack lengths were tested 
in order to assess the uniqueness of G c , the critical strain energy 
release rate, as the initial crack length was varied. Figs. 3.11 and 
3.12 show the impact response of a specimen with a 0.5 inch initial 
delamination, and the measurements taken of the resulting crack 
extension. Similar data for a specimen with a 2.0 inch initial 
delamination are shown in Figs. 3.13 and 3.14. In both cases, the 
impact velocity is very close to the threshold velocity required to 
cause crack extension in specimens of that particular geometry. 
Accordingly, significant variations in the threshold velocity from 
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that determined for the configuration with a 1.0 inch long initial 
delamination (Fig. 3.3) are noted. It would be expected, based on 
elementary fracture mechanics, that the impact velocity required to 
initiate crack propagation in a laminate would decrease as the length 
of the initial crack is increased. This trend is reflected in the 
values of threshold velocity for the specimens with 0.5, 1.0, and 2.0 
inch initial delaminations. The effect of initial crack length on G c 
is yet to be established, however. 

A comparison of the dynamic strain energy release rate prior to 
crack extension in specimens with 0.5, 1.0, and 2.0 inch initial 
delaminations should indicate if is independent of the initial crack 
length. G^ may be characteristic of the material, analogous to the 
fracture toughness in static loading. An analysis of these three 
cases is presented in Chapter 5. 

3.2.2 Off-Midplane Delamination 

All of the cases discussed so far involved delamination along the 
midplane of the beam. If the embedded crack is placed at a different 
through-the-thickness location, different crack propagation character- 
istics are observed. In the following impact specimens, the embedded 
crack is halfway between the beam midplane and outer surface. Thus, 
five plies are on one side of the crack and fifteen on the other. For 
these specimens, the camera was oriented to record the propagation of 
both crack tips simultaneously, instead of only a single crack tip, as 
in the previous cases. 

Some distinctly different features of crack propagation in this 
case can be seen from the photographs in Figs. 3.15 and 3.17. 
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"Buckling" of the delaminated plies is seen to occur at 125, 813 and 
875 ys in Fig. 3.15 and at 63, 813 and 875 ys in Fig. 3.17. The 
* photographs suggest, then, that the primary mode of crack extension in 
this case involves more Mode I (opening) than Mode II (shearing) type 
of action. The delamination buckling phenomenon has been noted in 
numerous other applications involving static, dynamic, and fatigue 
loading [28-35]. 

Tables 3.2 and 3.3 show that considerably greater impact energy 
is required to initiate crack propagation when the embedded crack lies 
along the midplane. The fact that no crack opening similar to that 
shown for off-midplane cracks is seen for midplane cracks (Figs. 3.7 
and 3.9) suggests that considerably less Mode I action is involved when 
the crack lies on the midplane. 

The intermittent nature of the delamination process is illustrated 
in Figs. 3.16 and 3.18 after the onset of crack propagation has 
occurred. Flexural wave propagation through the delaminated plies 
causes them to exhibit a beam-like dynamic behavior independent of 
the gross deformation of the specimen. Reflection of the waves between 
crack tips causes alternating propagation-arrest of the crack tips 
similar to that shown in Fig. 3.18 and to a lesser extent in Fig. 3.16. 

3.3 Stiffness Loss Due to Delamination 

Because the natural frequencies of a structure are dependent on 
its stiffness, a decrease in stiffness such as that caused by delamina- 
tion should cause a corresponding decrease in the natural frequencies. 
To investigate this effect, the first five natural frequencies of a 
series of midplane-cracked specimens with geometry similar to 
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Configuration B in Fig. 3.1 were determined with the apparatus shown 
in Fig. 3.19 before and after they were impacted. 

A sinusoidally time varying force was applied to the free end of 
the cantilevered specimens by means of a small magnet of negligible 
mass which was glued to the free end of the specimen, and a stationary 
magnetic exciter connected to a wave form generator. A Hewlett-Packard 
HP-141T Spectrum Analyzer slowly varied the forcing frequency of the 
wave form generator. This created a sinusoidally time varying magnetic 
field that was used to apply the harmonic force to the specimen 
through the magnet attached to its free end. A microphone was used to 
measure the magnitude of the output response (the motion of the 
specimen) accoustical ly. This was displayed graphically on the 
Spectrum Analyzer versus the input frequency. A dual trace oscilli- 
scope was used to display the frequency and magnitude of both the 
input and output graphically. 

The ultrasonic C-Scans of several typical impacted specimens 
shown in Fig. 2.5 indicate that the variation of crack length through 
the width of the specimen is small. Therefore, a two dimensional 
finite element analysis, which will necessarily assume a uniform 
through-the-width crack, should accurately predict the vibration 
characteristics of the damaged specimens. Four-node isoparametric 
plane-strain finite elements [36] were used to model the damaged and 
undamaged specimens. Figures 3.20-3.24 compare the measured and 
calculated results for the variation in the first five bending 
frequencies as a function of crack area. These results indicate that 
two dimensional elements can accurately model the low vibration modes 
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of both damaged and undamaged specimens of this geometry 
the results seem to show that the higher frequencies are 
more sensitive to delamination damage. 


In addition, 
progressively 
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TABLE 3.1 Variation of Delaminated Area with Impact 
Velocity for Specimen Configuration A 
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*initial delaminated area is 1.00 in2 







TABLE 3.2 Variation of Delaminated Area wi 
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TABLE 3.3 Variation of Delaminated Area with Impact 
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TABLE 3.5 Variation of Delaminated Area with Impact 
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TABLE 3.7 Variation of Delaminated Area with Impact 
Velocity for Specimen Configuration H 
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Figure 3.1 Impact Specimen Configurations 
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Figure 3.2 Delaminated Area Versus Impact Energy 
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Figure 3.3 Delaminated Area Versus Impact Energy 
for Specimen Configuration B 




Figure 3.5 Delaminated Area Versus Impact Energy 
for Specimen Configuration D 
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Figure 3.7 Delaminated Area Versus Impact Energy 
for Specimen Configuration H 
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Figure 3.20 Vibration Test Set-Up 
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4. MODELING THE IMPACT FORCE 


In order to perform a finite element analysis of the impacted 
composite specimen, the dynamic force history due to the impact must 
be determined. In lieu of a direct measurement of the contact force 
between the impactor and the target specimen, several methods were 

used to estimate the actual force histbry. These will be described 
here. 

4.1 Model 1 - Preliminary Model 

Daniel et al. [37] conducted an impact experiment on boron/epoxy 
and graphite/epoxy composite laminates using a 0.3125 in. diameter 
silicon rubber ball as impactor. Although the contact force was not 
measured, they were able to determine the contact 'area as a function 
of time. The contact area versus time curve was approximated by a 
sine function. Although the exact relation between the contact force 
and contact area is still unknown, it seems reasonable to assume that 
the contact rorce can also be approximated by a sine curve as 

F(t) = F 0 sin (^) 0 < t < T 

=0 t > T (4 * 1} 

where T is the contact duration. To determine the unknown coefficients 
F q and T following experiment was performed. 
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An uncracked cantilever beam specimen, shown schematically in 
Fig. 4.1, was impacted with the half-inch diameter silicon rubber ball 
at the velocity of 1166 in/sec. Two strain gages (Micro Measurements 
ED-DY-031CF-350, S g * 3.25) were mounted on the back side of the 
specimen to measure the bending strain history. One of the gages was 
mounted directly opposite the impact point, and the other gage was 
placed at 2.0 inches away from the first gage. The strain histories 
measured by these two gages are presented in Figs. 4.1 and 4.2. 

A finite element model using 4-node isoparametric plane strain 

finite elements was then used to model the impacted beam and calculate 

the strains at the two gage locations. A uniform mesh of 1200 elements 

(2800 d.o.f.) was found to yield a converged solution and was used to 

find the values of F Q and T that best matched the experimental results. 

A detailed description of the finite element modeling procedure used 

is given in the next chapter. The finite element results shown by the 

dashed line in Figs. 4.1 and 4.2 were obtained with F Q = 44 lb. and 

T * 225 ys. In fitting these values, it was found more convenient 

to vary T to fit the time-phase and then determine the force amplitude 

F , as the strain is linearly proportional to the amplitude, 
o 

4.2 Model 2 - Experimental Model 

The dynamic response of the impact specimen will be predicted 
more accurately if the estimate of the force history is improved. The 
actual force history from the impact of the rubber impactor was 
therefore determined experimentally. 
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4.2.1 Impact of a Longitudinal Bar 

Longitudinal stress waves propagate nondispersively in a uniform 

thin bar. An instrumented bar is therefore well suited to determining 

the force history due to longitudinal impact [10, 38]. It is 

conceivable, if the deflection of the target does not significantly 

« 

affect the contact force, that the force history determined from the 
impact of the longitudinal bar would be nearly the same as that 
experienced by the composite laminate under the same impact conditions. 

The experimental set-up used is shown schematically in Fig. 4.3. 
The equipment used is similar to that used for the impact of the 
cracked laminates with the addition of the apparatus needed for 
measuring the strain history. As the strain pulse passes the gage, 
the measured change in voltage output is amplified by the pre-amp. 

The resulting strain history is stored in the Biomation waveform 
recorder and displayed on the oscilloscope. If necessary, the strain 
history can be plotted on paper using the X-Y plotter. The force 
history can be calculated directly from the measured strain. 

4.2.1. 1 Hertzian Contact Law - Analysis of Low Speed Longitudinal 
Impact of a Bar 

Hertz [39, 40] derived the force-indentation relation for the 
general case of two spherical elastic masses coming in contact with 
each other. To summarize the Hertzian contact law, we have 

F = Ka n (4.2) 

where 
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F = contact force between spheres 
a = relative indentation between spheres ( v i _v 2 ) 
n = 1.5 


and 



(4.3) 


where 


R.. = radii of spheres 



and Ej, are the respective elastic constants. A special case that 
is of interest here occurs when the target is flat (f^ * ®) in which 
case eq. (4.3) simplifies to 



Before proceeding to the more complicated problem of the high 
speed impact of the rubber ball, a preliminary analysis of the low 
speed impact of a steel ball on the aluminum bar was performed. A 
finite element program was developed to analyze the longitudinal 
impact of a bar. Four degree-of-freedom truss elements [41] were used 
to model the bar. 

The equations of motion for a structure subjected to a time 
dependent force F(t) are 
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{ F( t) } = [K] {6} + [M] {6} 


(4.5) 


' where [K] is the stiffness matrix, 

[M] is the mass matrix, 

{6} is the displacement vector, and 

{6} is the acceleration vector. 

In the formulation of this element, both the displacements and the 
strains at each nodal point are used as degrees of freedom. For the 
truss element, then, we have 


and 


and 


{6}^ = [u r (au/axjj, u 2 , (3 u/3x) 2 ] 


(4.6) 


[K] = EA 


6/5L 



sym. 

1/10 

2L/15 


. 

-6/5L 

-1/10 

6/5L 


1/10 

-L/30 

-1/10 

2L/15 


(4.7) 


[M] = 


AL 


420 


156 



sym 

22L 

4L 2 
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13L 

156 


-13L 

-3L 2 

-22L 

4L 2 
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where 


are the nodal displacements, 

EA is the axial stiffness of the element, 
pA is the mass per unit length, and 
L is the length of the element 

A detailed derivation of the equations of motion for this high- 
order truss element is given in [41]. 

The only non-zero term in the force vector ( F( t ) } is that 
corresponding to the displacement at the impact point on the bar. The 
contact behavior between relatively hard materials such as steel and 
aluminum is well described by the Hertzian law. The Hertzian contact 
law was therefore incorporated into the longitudinal bar finite 
element program to perform the low speed longitudinal impact analysis. 

The Hertzian impact force is given by eq. (4.2). In this case 
the indentation is 


o * x-u 


(4.8) 


where 


x = displacement of the impactor 
u = bar displacement at the impact point 

The equation of motion for the center of mass of the impactor is 


F = mx 


(4.9) 


where F is the impact force given in eq. (4.2) and m is the mass of the 
impactor. 
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Equation (4.9) must be solved simultaneously with the equations 
of motion for the bar to determine the dynamic force history. The 
* initial conditions are given as 


x ( o ) = x(o) = u(o) = F(o) = 0, x(o) = v 

(4.10) 

where v is the impact velocity. The following algorithm was 
determine the impact force. 

used to 

x m * x i + x i 4t 4 x i H 1 - 

(4.11) 

F i+1 * k * x i+l " u j>" 

(4.12) 

x 1+l * - F i + l' M 

(4.13) 

x i+l ‘ *1 + X 1 4t 

(4.14) 


where At is the integration time step and the notation x i * x(t = iAt) 
has been used. Equation (4.12) is only an estimate of F. + j because the 
displacement u^ + j at the impact point has been estimated by from the 
previous time step. 

The experimental set up for the low speed impact tests is 
identical to that shown in Fig. 4.3 with the exception that the steel 
impactor was suspended as a pendulum and dropped, instead of being 
fired ballistically with the air gun. Figure 4.4 compares the measured 
force history with that obtained using a 39-element (80 d.o.f.) model 
of the bar and an integration time step At of 0.5 ysec in the finite 
element program. Figure 4.5 shows the strain history measured at a 
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single gage location on the bar compared to the finite element result. 
These results seen to verify that the longitudinal bar finite element 
•program accurately predicts the strain history due to a Hertzian 
impact. Figure 4.6 shows the propagation and reflection of the strain 
pulse along the bar as calculated by the finite element program. The 
boundary conditions of the bar are free-free, so the compressive 
incident pulse reflects from the distal end of the bar as a tensile 
pulse of identical shape. 

4. 2. 1.2 High Speed Longitudinal Impact of a Bar 

The ballistic impact set up shown in Fig. 4.3 was used to measure 
the strain history resulting from impact of the rubber ball on the 
aluminum bar. The force history can be calculated directly from the 
measured strain. 

Figure 4.7 shows a typical measured impact force versus time 
behavior for the impact of the half-inch diameter ball. A simple 
approximation can be used to describe the curve: 


* 


F(t) = 


fFo sin<#-) 
r o 

*( t-tp ) 

^ F o cos 2(T-tp°J 
o 


O < t < tp 

r o 


tp < t < T 
r o 

t > T 


(4.15) 


where F is the maximum force, tp is the time when the maximum force 
occurs, and T is the contact duration. Each of these three parameters 
used in describing the approximate force history can be read directly 
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from a three-parameter estimate of the actual force history such as 
that shown in Fig. 4.7. 

A series of tests was performed to determine how the force history 
varied with impact velocity. Figures 4.8-4.10 show the variation in 
shape and amplitude of the force history resulting from impact of the 
half-inch ball on the bar. Figure 4.8 shows that the amplitude of 
the force varies in proportion to v^, which is in contrast to the 
linear behavior predicted by the simple spring-mass single-dof impact 
model [42]. Figure 4.9 shows that the contact time varies inversely 
with the impact velocity, which, is also in contrast to the spring-mass 
model, which predicts that the two are independent. Figure 4.11 shows 
that the impulse measured from the experimental data varies linearly 
with impact velocity. 

A second series of tests was conducted with a smaller (3/8 inch 

diameter) ball of the same material in order to further characterize 

the impactor. The results presented in Figs. 4.12-4.15 show that the 

trends followed by the data are similar to those observed for the 

impact of the larger ball. The maximum force varies in proportion 
2 

to v , the contact duration varies inversely with v, and the shape 

tfr </"T of the force history is approximately constant over the velocity 
range tested. 

4.2.2 High Speed Transverse Impact of a Laminated Beam 

The impact force history for the specimen shown in Fig. 4.1 is 

adequately approximated by eq. (4.15) with the parameters F = 42 lb 

o * * 

tp - 35 psec, and T = 415 psec. The calculated strain response is 
compared with the measured values in Figs. 4.1 and 4.2. There is a 
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significant improvement over the results obtained with the two- 
parameter approximation in eq. (4.1). Therefore, in all subsequent 
-impact analyses, the force history used will be of the form given in 
eq. (4.5). 

To establish Force History versus Impact Velocity relations for 
the beam specimen analogous to those shown in Figs. 4.8-4.10 for the 
bar, a further series of experiments was conducted with composite 
laminates of geometry and lay-up similar to those used in the crack 
initiation studies. The notable differences in these specimens were 
that they had no initial delamination, and two strain gages 
(MM ED-DY-031CF-350, S g = 3.25) were mounted on each of them as shown 
schematically in Fig. 4.16. Two specimens of different lengths were 
tested in order to see the effect of the specimen size on the force 
history. The strain history at both gage locations was recorded up 
to approximately 1000 ys after impact with the equipment shown in Fig. 
4.3. Strain histories were obtained for impact velocities ranging 
from 1000 to 4000 inches per second. 

The composite laminate was modeled with plane strain finite 
elements. A uniform mesh of 1200 elements and 2800 total degrees of 
freedom was found to yield a converged solution which adequately 
modeled the dynamic response of the laminate. A more detailed 
description of the finite element modeling procedure is given in 
Chapter 5. The three parameters in the assumed force history given 
in eq. (4.5) were varied in the finite element analysis to match the 
measured strain response. Comparisons of the measured and numerically 
obtained strains are shown in Figs. 4.16-4.23 for the shorter specimen 
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and in Figs. 4.24-4.31 for the longer specimen. The different 

parameters used in the force histories for these analyses are given in 

Tables 4.1 and 4.2. The variation of each of the three parameters 

with impact velocity is shown graphically in Figs. 4.32-4.34. Figures 

4.32 and 4.33 indicate that the contact force is lower and the total 

contact time is longer in the shorter specimen, for a given impact 

. velocity. This indicates that the flexural wave reflections from the 

boundaries have a significnat effect on the force history. If the 

effect of the wave motion was not considered, and a single degree of 

freedom spring-mass model based on the static stiffness of the 

respective laminates was used to anticipate these trends in the force 

** history, the opposite behavior would be predicted. Indeed, the shape 

of the force history varies with the length of the target as shown in 

Fig. 4.34. Comparing these results with those obtained from the 

longitudinal bar experiment shown in Figs. 4.8-4.10, it is apparent 

that the relationships between the impact velocity and the three 

parameters describing the force history are similar. That is, the 

maximum force varies in proportion to V 2 , the contact duration varies 

inversely with V, and the ratio tp /T is relatively constant over the 

o 

velocity range tested. Although both sets of data follow similar 
trends, it is also apparent that the different dynamic response of the 
targets significantly affects the variation of the force history with 
impact velocity. For this reason, the contact force measured in the 
bar experiment cannot be applied directly to the impact analysis of 
the laminates. For laminates of the specific dimensions tested, the 

data presented in Figs. 4.32-4.34 give more accurate estimates of the 
actual impact force history. 
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4.3 Contact Area and Impact Force 

It is evident from Fig. 3.7 that a large amount of deformation 
occurs in the rubber impactor while it is in contact with the 
composite laminate. This deformation spreads the load due to contact 
over a larger area than would occur with a more rigid impactor. In 
calculating the loading used in the finite element model, this force 
distribution must be accounted for if accurate predictions of the 
strain near the impact point are needed. 

For the purpose of modeling the variation of contact area with 
time, it is assumed that 

MU = (4.16) 

A o F o 

where A(t), F(t) are contact area and contact force, respectively, and 

A , F are their maximum values. The actual circular contact area of 
o o 

radius r is approximated by a rectangular strip of dimension 2r x 1. 
Thus, 

r(t) " r °JW (4 ' 17> 

where r(t), r Q are the contact length and its maximum value. The 
maximum contact area can be measured from the imprint left by the 
impactor on the laminate. It is further assumed that the spatial force 
distribution f(r) over the contact length at any given time is as shown 
in Fig. 4.35. Hence, 
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(4.18) 


f(r) = F max cos <f r < r 0 


The specific nodal force distribution is shown schematically i 
4.35 and is written as 


r. + Ar/2 

f i = f T(r) dr 

} r. - Ar/2 


After imposing the condition 


+r 

f nr) 

' -r 


dr = F 


where F is the total force at the given time, (4.5) becomes 


f i = I tsin 4r^ {2r i + Ar ) - sin ^ (2r. - Ar)] 


n Fig. 


(4.19) 


(4.20) 


(4.21) 
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TABLE 4.1 Force Histories Used in Analysis of Impact 
on Short Composite Beam Specimen 


Force History 

V(in/s) 

- .. . - -J 

F„Ob) 

t F (ys) 
0 

T(ps) 

r 0 On) 

1 

1166 

42 

35 

415 

0.25 

2 

1981 

84 

20 

310 

0.25 

3 

2994 

170 

13 

205 

0.30 

4 

4049 

245 

10 

155 

0.35 










TABLE 4.2 Force Histories Used in Analysis of Impact 
on Long Composite Beam Specimen 


Force History 

V( in/s) 

F o< lb > 

t p (ps) 
0 

T(ps) 

r „(in) 

1 

1166 

44 

30 

355 


2 

1932 

90 

21 

250 

■ 

3 

3058 

' 215 

. 13 

I 158 


4 

4026 

325 

10 

120 

0.35 
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Figure 4.1 Strain History at Impact Point of Beam Subjected 
to 1166 in/s Impact with Half-Inch Diameter 
Impactor 
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Figure 4.2 Strain History 2.0 Inches From Impact Point 
of Beam Subjected to 1166 in/s Impact with 
Half- Inch Diameter Rubber Impactor 
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Figure 4.4 Force History From 62. lin/s Impact of 

5/8 Inch Diameter Steel Ball on Aluminum 
Bar 
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Figure 4.5 Strain History at Midpoint of Bar From 
62.1 .in/s Impact of 5/8 Inch Diameter 
Steel Ball 
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Figure 4.6 Propagation and Reflection of Longitudinal Strain 
Pulse From 62.1 in/s Impact of 5/8 Inch Diameter 
Steel Ball 
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Figure 4.7 Force History Measured From 3500 in/s Impact 
of Half-Inch Diameter Impactor 
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Figure 4.10 t R /T Versus Impact Velocity for 
o 

Half-Inch Diameter Impactor 
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Figure 4.11 Measured Impulse Versus Incident Momentum 
for Ha If- Inch Diameter Impactor 
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o 

Figure 4.12 Force Amplitude Versus V for 3/8 Inch 
Diameter Impactor 
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Figure 4.14 t p /T Versus Impact Velocity for 3/8 Inch 
♦ o 

Diameter Impactor 
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Figure 4.15 Measured Impulse Versus Incident Momentum for 
3/8 Inch Diameter Impactor 
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Figure 4.17 Strain History at Gage 2 in Short Beam 
Specimen Using Force History No. 1. 
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Figure 4.20 Strain History at Gage 1 in Short Beam 
Specimen Using Force History No. 3. 
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Figure 4.21 Strain History at Gage 2 in Short Beam 
Specimen Using Force History No. 3. 
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Figure 4.22 Strain History at Gage 1 in Short Beam 
Specimen Using Force History No. 4. 
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Figure 4.23 
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Strain History at Gage 2 in Short Beam 
Specimen Using Force History No. 4. 
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Figure 4.24 Strain History at Gage 1 in Long Beam 
Specimen Using Force History No. 1. 
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Figure 4.25 Strain History at Gage 1 in Long Beam 
Specimen Using Force History No. 1. 
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Figure 4.26 Strain History at Gage 1 in Long Beam 
Specimen Using Force History No. 2. 
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Figure 4.27 Strain History at Gage 2 in Long Beam 
Specimen Using Force History No. 2. 
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Figure 4.28 Strain History at Gage 1 in Long Beam 
Specimen Using Force History No. 3. 
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Figure 4.29 Strain History at Gage 2 in Long Beam 
Specimen Using Force History No. 3. 
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Figure 4.31 Strain History at Gage 2 .in Long Beam 
Specimen Using Force History No. 4. 
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Figure 4.33 Inverse Contact Time Versus Impact Velocity 
for Impact on Composite Laminates. 
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Figure 4.34 tp /T Versus Impact Velocity for Impact 
on°Composite Laminates. 
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Figure 4.35 Nodal Force Distribution at a 
Fixed Instant in Time 
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5. INSTABILITY OF DELAMINATION CRACKS - 
CRITICAL STRAIN ENERGY RELEASE RATE 


Strictly speaking, the impact problem we are concerned with is a 
three-dimensional problem. However, photographs taken by the high 
speed camera indicate that the impactor deformation covered almost the 
whole width of the composite beam specimen. Moreover, due to the small 
dimension in width, the specimen behaved like a beam except during the 
initial period of contact. In view of the foregoing, the laminate 
specimen was approximated as a two-dimensional body and a two- 
dimensional linear elastic finite program was used to perform the 
dynamic analysis. The impact load was taken to be uniform across the 
width of the specimen, and a state of plane strain parallel to the 
longitudinal cross-section was assumed. This cross-section was then 
modeled by regular four-node quadrilateral isoparametric finite 
elements. 

5.1 Equivalent Moduli 

Ideally, each lamina should be modeled with a number of finite 
elements to ensure the best accuracy. However, such a procedure would 
lead to a formidably large number of elements for the 20-plied laminate. 
For this reason, the [90/0]g s laminate was transformed into an equiva- 
lent homogeneous plate with a set of effective moduli obtained by 
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using appropriate constant strain and constant stress assumptions 
[43]. The equivalent material properties derived below are used to 
model the laminate in regions remote from the crack tip. 

The material properties for the individual orthotropic lamina are 
defined relative to the coordinate system shown in Fig. 5.1 and are 
defined as 


ij 


1 3 


Young's Modulus in the * i 1 direction 
i - 1, 2, 3 

Shear Modulus in the i-j plane 
i, j - 1, 2, 3 G.j = G j i 

-ej/e.j for uniaxial stress = o 
all other stresses are zero 
i, j = 1. 2, 3 and 


lii = !ii 

E i 


(5.1) 


The equivalent material properties for the [0/90] laminate with 
reference axes defined as shown in Fig. 5.2 will be given the 
corresponding notation E.^ G-, and v.. and are derived as follows. 

I J 1 J I J 

Assume the laminate in Fig. 5.2 is subjected to a uniaxial stress 
Oj^ in the 1-direction, and that it will undergo a resultant uniform 
extension in the two plies 


90 0 

e ll = e ll * e ll 


(5.2) 


where e^ is the engineering strain in the * i ' direction of the 
laminate and e^ is the engineering strain in the * i ' direction of the 
e-oriented ply. 
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Under this condition, the stress is not the same in both 
The effective stress for the laminate is given by 


11 


90 .90 . 0 0 

. °n A * °n A 

A 90 + A° 


where 


°i -j is the laminate stress in the * i 1 direction 
0 

a i -j is the stress in the e-oriented ply in the 
laminate 'i 1 direction 

0 

A is the cross sectional area of the 6-oriented 

ply. 


Using the volume fractions 

A 90 A 0 

V 90 * ft 90 ~ A 0 V o * “90 + A 0 

eq. (5.3) can be written 

°11 = °11 V 90 + a ll V 0 
The effective Young's Modulus is defined by 

! i = o ii /e n 

since 


plies. 


(5.3) 


(5.4) 


(5.5) 


(5.6) 
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(5.7) 


_ r p a 

°11 " t 2 e ll t 2 e ll 

and 

0 _ r = Ee 
°11 t l e ll t l e ll 

Using eqs. (5.7) and (5.8), eq. (5.5) becomes 

°11 = ^ E 2 V 90 + E 1 V 0> e ll 

so, from eq. (5.6) 


(5.8) 


(5.9) 


E 



90 


+ E, V 


1 v o 


(5.10) 


in this case, V^q = Vq = 1/2, so 


= E 1 + E 2 

4 2 


(5.11) 


for the [0/90] laminate. A similar procedure can be used to show 


t 2 c r 

Assume now that the laminate in Fig. 5.2 is subjected to a 
constant stress in the 3-direction. Assume 


90 _ 0 

°33 ' °33 ' 33 


( 5 . 12 ) 
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Now 


90 /c 

e 33 °33 /E 3 


(5.13) 


and 


e 33 = °33 /E 3 


(5.14) 


The strain in the laminate is given by 


> = e 33 t 90 * e 33*0 

33 *90 + *o 


(5.15) 


where t Q is the thickness of the e-oriented ply. Analogous to eq. 
(5.4), the volume fractions can be defined as 


90 


90 *90 + *0 


v = y 

0 *90 + *0 


(5.16) 


so eq. (5.15) can be written 


*33 = e 33 V 90 + e 33 V 0 

= I fe 90 + e° 1 
2 le 33 e 33' 


(5.17) 


Using eqs. (5.13) and (5.14), eq. (5.17) becomes 


e 33 = °33 /E 3 = °33 /E 3 


(5.18) 


so 
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(5.19) 


= E 


3 


0 

for the [0/90] laminate. 

Assume now that the laminate in Fig. 5.2 is subjected to a 
shearing stress t 13 on the 1-face and oriented in the 3-direction. 
The shear stress must be continuous at the ply interface, so 


90 _ 0 

t 13 ' t 13 ‘ t 13 


The shear strains are 


90 _ Il3 v ° 111 

Y 13 ‘ G 23 Y 13 S 13 


(5.20) 


(5.21) 


Defining as the shear strain of the laminate. 


r 13 


90 t + v 0 t 

Y 13 *90 y 13 t 0 

*90 + *0 


(5.22) 


Using eqs. (5.16), eq. (5.22) becomes 


90 

y 13 “ y 13 


'90 


y 13 V 0 


1 , 90 

2 (y 13 


* n°3> 


(5.23) 


Using eqs. (5.21) and assuming 

G 23 * G 13 


(5.24) 


eq. (5.23) becomes 
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y 13 " t 13 /G 13 = t 13 /G 13 


(5.25) 


so 


G 13 ( = G 13 


for the [0/90] laminate. 

Assume now that the laminate in Fig. 5.2 undergoes 
strain e^ in the 1-direction. Assume 


90 0 

11 " e ll “ e ll 


The strain in the 90-degree ply is given by 


.90 


'33 ' “ v 23 e ll 


and in the 0-degree ply is 


e 33 * “ v 13 e ll 


The laminate strain is 


'33 


e 90 t + e° t 
_ e 33 t 90 e 33 1 0 


*90 + *0 


Using eq. (5.16), eq. (5.30) becomes 


(5.26) 


uniform 


(5.27) 


(5.28) 


(5.29) 


(5.30) 
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'33 


e 90 w + gO v 
e 33 90 33 


- fe 90 + e u ) 

2 * e 33 e 33 ; 

\ ( v 9i + v j3) e n 


0 


'23 


(5.31) 


Now 


e 33 = ~ v 13 e ll 


(5.32) 


so 

v 13 = \ ^ v 13 + v 23^ 


(5.33) 


for the [0/90] laminate 

The mechanical properties of the T300/934 graphite epoxy lamina 
are given as 


Ej = 19.5 x 10 6 psi 

E 2 * 1.5 x 10 6 psi 

G 12 - 0.725 x 10 6 psi 

v 12 = v 13 “ v 23 s 0,33 
In addition it was assumed that 

G 23 " G 13 = G 12 
and E 2 = ^2 


(5.34) 


(5.35) 
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in the numerical analysis. Corresponding effective Moduli derived 
previously are 

Ej = = 10.5 x 10^ psi 

E^ = 1.5 x 10^ psi 

G 13 = 0.725 x 10 6 psi (5.36) 

~ v 23 = 0*33 9 v j2 = *025 
P = 1.49 x 10" 4 lb-s 2 /in 4 

These equivalent moduli are used as the elastic constants of the 
finite elements in regions remote from the crack tip. 

5.2 Finite Element Modeling 

Four-node quadrilateral isoparametric elements [44] were used to 
model the composite laminates. A state of plane strain is assumed in 
the y-z plane so the displacement field is given by 

v = v(y , z) 

w = w (*» z ) (5.37) 

In the isoparametric formulation, the element shape functions are 
used to express the displacement field within the element in terms of 
the nodal displacements 

4 

v(y, z) - £ N. v. (5.38) 


and 
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4 

w(y, z) = l N. w. 
i=l 1 

where v., w i are the nodal displacements, and N i are the element shape 
functions, given by 

Ni = \ (1 - 00- n) 

N« = j (i + 0(1 - n) 

(5.39) 

N 3 = \ (1 + 00 + n) 

N 4 = \ (1 - 0(1 + n) 

and c, n are the curvilinear coordinates used to define the element 
boundaries . 

Using the C, n coordinate system, the rectangular coordinates x,y 
are written in terms of the nodal coordinates as 

4 

y - Jj N i 

and < 5 ' 40 > 

2 * I "i 

i=l 1 1 

So the coordinates of nodes 1,2,3, and 4 of the quadrilateral element 
are (-1, -1), (1, -1). (1. D ("1» D respectively, when written in 
terms of the curvilinear coordinates (c, n). 
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The effective elastic constants calculated in Section 5.1 were 
used as the material properties of the finite elements in regions 
remote from the crack tip. The through-thickness distribution of 
elements used in the remote regions is shown in Fig. 5.3. The exact 
material properties given in eqs. (5.34) were used to model the 
laminate in the region near the crack tip, where an accurate represen- 
tation of the steep stress gradient is crucial. Figure 5.4 shows the 
transition from the coarse mesh in the remote regions to the finer 
mesh near the crack tip. The elements in the central 0.1365 inch long 
region of the mesh were modeled with the exact single-ply material 
properties given in eqs. (5.34). The elements outside of this region 

were assumed to have the effective material properties derived in the 
previous section. 

An accurate finite element representation of the regions remote 
from the crack tip is necessary to adequately model the motion of the 
laminate. Therefore, an investigation of the effect that the size of 
the remote elements has on the calculated response at the crack tip 
was conducted. For this study, all elements were assumed to have the 
equivalent material properties calculated earlier. Thus, the only 
effect on the crack tip energy calculation was that due to mesh 
refinement, not variability in the assumed material representation. 

In addition, the near tip mesh was relatively coarse, so the system 
of equations would not become prohibitively large as the remote mesh 
was refined. For this study, the near tip region was modeled with six 
elements through the beam thickness, as shown in Fig. 5.3, and an 
assumed virtual crack extension ratio aa/a of 1/100. Figure 5.5 shows 
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the effect of remote mesh refinement on the calculated energy release 
rate. If the mesh is too coarse, the finite element model is 
* excessively stiff, and the energy release rate is over estimated. 

When the element aspect ratio is reduced to 4:1, the solution has 
apparently converged. This is indicated by the insignificant change 
in the calculated solutions for the 4:1 and the finer 2:1 aspect ratios. 
Therefore, an aspect ratio of 4:1 is considered the maximum allowable 
value required to accurately model the bending of the beam in regions 
remote from the crack. 

A second convergence study was performed for the mesh near the 
crack tip. If this mesh is too coarse* the stress gradient near the 
crack tip will not be adequately represented, and the strain energy 
release rate calculation will be affected. In this study, the aspect 
ratio of the remote mesh was kept at 4:1 while successive refinements 
of the near-tip mesh were performed. Figure 5.6 indicates that the 
finite element models using 20 and 40 elements through the beam 
thickness represent the near-tip response with equal accuracy. A 
small but consistent variation from the solution obtained using the 
coarsest mesh is evident. The two convergence studies have thus 
indicated that a finite element mesh with a maximum 4:1 aspect ratio 
remote from the crack and 20 elements through the beam thickness in 
the region nearest the crack tip will yield a converged solution for 
the beam with homogeneous material properties. This mesh convergence 
will also be assumed valid for the laminated beam. 

The quadrilaterial elements were taken to be square near the crack 
tip. Thus, for the twenty-plied laminate the virtual crack extension 
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Aa used in the energy release rate calculation and assumed equal to 
one near-tip element length [25], is equal to the thickness of one 
lamina, approximately .005 inch. 

Newmark's implicit time integration scheme [45] was used to 
numerically integrate the equations of motion. For the choice of 
parameters used in the integration, y = 0.5 and 6 = 0.25, this method 
is also known as the "Constant Average Acceleration" method, and is 
unconditionally stable for linear problems. Thus, a relatively large 
time step. At = 2.5 usee, could be used to obtain a converged solution. 

The three-parameter model of the impact force history, discussed 
in the previous chapter, was used in the following analyses of the 
cracked laminates. The particular parameters used for these cases 
were extrapolated from the data in Figs. 4.32-4.34. The force 
histories for each of the three following analyses are given in 
Table 5.1. 

5.3 Strain Energy Release Rate 

Of interest to the present study is finding a parameter that can 
be used to gage the on-set of dynamic delamination crack propagation. 

A natural choice is the use of dynamic strain energy release rate 6, 
which can be calculated by using the crack-closure energy given by 
[24] 


1 f Aa 

G ' la?o 255 Jo V * V )dx < 5 - 41 > 

in which o yy and c xy are evaluated at the original crack size a, and 
u and v are the near-tip displacements corresponding to the extended 
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crack of length (a + Aa). Using the finite element method, the 
integral in eq. (5.41) can be carried out by using discrete nodal 
' forces and displacements. Moreover, if a fine mesh is used, i.e. 

Aa « a, then crack opening displacements u and v can be approximated 
by those for crack of length. The finite element mesh near the crack 
tip of interest is shown schematically in Fig. 5.7. Using the crack 
closure method of calculating the strain energy release rate [25], the 
Mode I and Mode II contributions can be calculated separately: 

G i = ’ 1ra Si f a<V v b> 

Aa+o 

and (5-«) 

G n • si Vv u b> 

Aa-*-o 

where T and F are the nodal forces in the u and v directions that are 
d 9 

needed to "close" the crack by an amount Aa. For a small enough 
element length Aa, these can be approximated by T c and F c , respectively. 

5.3.1 Verification of Crack Closure Method 

A centrally cracked rectangular panel of homogeneous isotropic 
material subjected to a uniform tensile step function loading was 
analyzed by Chen using a finite difference method [46]. His solution 
was used in this study to validate the aforementioned finite element 
method in conjunction with the crack closure energy calculation. To 
compare with Chen's solution, which was presented in terms of stress 
intensity factors, the following relation for Mode I fracture 
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V‘> - fr K i 2(t) 


(5.43) 


was used. In eq. (5.43), y is the shear modulus, Kj is the Mode I 
stress intensity factor and 

{ 3 - 4v for plane strain 

(3-v)/(l+v) for plane stress. (5.44) 

This relation can also be applied to stationary cracks under dynamic 
leading [47]. 

Figure 5.8 shows the geometry and material constants of the model 
studied by Chen [46]. Due to symmetry, only a quadrant was modeled. 
Figure 5.9 shows the histories of the normalized stress intensity 
factor Kj, given by 



P/iTa 


(5.45) 


obtained by [46] and by the present method. 

Three finite element meshes were used. The coarse mesh consists 
of 99 4-node quadrilateral plane strain elements and 221 degrees of 
freedom. In the critical area near the crack tip, the mesh size 
yields a ratio of Aa/a = 1/4. The finer mesh is composed of 323 
elements with 682 degrees of freedom and a near-tip mesh size of 
Aa/a = 1/16. The result from the third mesh was found to agree very 
well with that from the second mesh and thus can be considered a 
converged solution. The comparison presented in Fig. 5.9 shows that 
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the present method is quite acceptable. 

5.4 Impact Analysis of Cracked Laminates 

The threshold impact velocities required to cause crack extension 
in each of the impact specimen configurations shown in Fig. 3.1 were 
determined from the data presented in Figs. 3. 2-3. 7. Finite element 
analyses were performed to calculate the dynamic strain energy release 
rate that results from an impact of each specimen configuration at its 
particular threshold velocity. The force history used in the analyses 
were determined from linear extrapolations of the data shown in Figs. 
4.32-4.34. The parameters obtained are given in Table 5.1. Because 
each analysis is performed at the threshold impact velocity for the 
particular specimen configuration, the maximum value of energy release 
rate reached after impact can be considered to be the critical value of 
G necessary to cause crack extension. THe maximum values calculated in 
the following analyses are therefore given in Table 5.1 as estimates 
of G c . 

Specimen configurations B, G, and H have initial crack lengths of 
1.0, 0.5, and 2.0 inches, respectively. The strain energy release rate 
was calculated at both crack tips for each of these specimens. The 
results indicate that the crack tip nearest the fixed end always 
reaches a higher maximum value of energy release rate, and can there- 
fore be considered the "critical crack tip. Experimental evidence that 
this is indeed the case for specimen configuration C, in which the 
initial delamination is not on the beam midplane, is presented in Figs. 
3.17 and 3.19. In both cases, the crack tip nearest the fixed end 
extends first. This crack tip must therefore reach the critical energy 
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release rate before the free-end crack tip. Two notable exceptions to 
this generality are configurations D and F. An earlier investigation 
*[48] showed that in these two cases, a "peeling" action occurs in 
which the crack tip nearest the free end extends to the free boundary 
before the other crack tip moves, causing separation of the specimen. 
The free-end crack is therefore the critical one in these cases. With 
the exception of the three cases mentioned earlier, the energy release 
rate was calculated only at the critical crack tip. 

The values of G c given in Table 5.1 for the different specimen 
configurations analyzed in Figs. 5.10-5.17 indicate that G £ for the 
Mode I cracks is considerably lower than that for Mode II. This trend 
has been clearly identified in experimental work using static loading 
[9], The empirical method used to find the impact force for this 
analysis may have predicted a higher force than actually occurred in 
the case of the higher impact velocities. The data in Fig. 4.2 begins 
to deviate slightly from the linear approximation as the impact 
velocity increases. Therefore, the energy release rate calculations 
for those specimens impacted at the higher velocities, most notably 
configuration G, may be less accurate than the other estimates. 

As expected from the experimental results, the Mode of crack 
extension varies as the location of the delamination through the 
thickness of the laminate is varied. The calculated energy release 
rates indicate that the initial crack extension for the off-midplane 
crack geometries involves primarily Mode I action, while the specimens 
with the initial delamination on the beam midplane undergo only Mode II 
crack extension. Although there is a small Mode II contribution for 
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the off-midplane specimens before the peak value is reached, it is 
insignificant in comparison to Mode I when G = G c » In the analysis of 
• the midplane cracked specimens, the Mode I contribution was always 
insignificant in comparison to Mode II. 

Static measurements of the critical strain energy release rate 

necessary to cause delamination crack extension in Graphite/Epoxy 

laminates have been performed by several investigators, with relatively 

consistent results. O'Brien [8] measured the strain energy release 

rate at the onset of delamination in tensile tests of [±30/±30/90/90] $ 

laminates. The critical value of G for a delamination at the -30/90 

2 

interface was found to be 0.78 in-lb. /in . This value was then used to 
predict the onset of delamination in laminates with different layups. 

The results indicated that G c may be independent of the ply orientations 
along the delamination interface. Wilkins [9] used Double Cantilever 
Beam and Cracked Lap Shear specimens to generate pure Mode I and Mode 
II cracks. Two different layups were used for both the Mode I and 
Mode II tests. For the Mode I case, the first ply configuration used 
was [O^/D/Oj^] where D represents the initial delamination produced 
by embedding Kapton film between the two center plies of the 24-ply 
unidirectional ly reinforced laminate. The other layup tested for the 
Mode I case was [0 2 /90/0g/D, 90] $ where D, 90 represents the embedded 
delamination adjacent to a 90 degree ply located on the midplane of the 
DCB specimen. In the first layup used, the delamination extends along 
a 0/0 interface and in the second layup, along a 0/90 interface. Two 
different ply orientations were also tested for the shear specimen. 

The first was [0/r45 4 /0/D/0/90/0 2 ], so ten plies are on one side of the 
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del ami nation and four on the other. In this case the delamination is 

embedded at a 0/0 interface. The second layup tested was [0/145^ + 

. 452 /O/D/ 9 O/O 2 ] so the embedded delamination lies along a 0/90 interface 

Critical energy release rate values were found to be 0.5 and 0.8 

in-lb. /in for the Mode I and Mode II cases, respectively. In addition 

the G c values for each case were nearly the same at the 0/90 and 0/0 

interfaces. This is further evidence that G c may be independent of the 

ply orientations along the delaminating interface. 

It should be expected that the magnitude of the strain energy 

release rate necessary to cause crack extension in the case of dynamic 

loading would be higher than that in the static case. This trend was 

found to be true for steels [14]. The dynamic values of G estimated 

c 

here show a similar relationship to the static values, so in that 
sense some confidence can be placed in these results. 
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Table 5.1 Force Histories Used in Analysis of Pre-Cracked 
Composite Beam Specimens 
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Finite Element Constants : 
Fj = E 2 = 10.5x10® psl 

F 3 c 1.5xlO b psi 
6 j 2 = 0.725x10® psl 


Lamina Properties : 
= 19.5x10® psl 

E 2 * 1.5x10® psi 

®12 = 0-725x10® psi 


v 12 


* v 13 = v 23 = °- 33 
1.49xl0 ‘ 4 -^4? 


v 13 * v 23 * °‘ 33 
v 12 * 0.05 


Figure 5.3 Finite Element Discretization of [90/0] c 
Graphite/Epo*y Laminate Remote from 5s 
Crack Tip 
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0.45 inch 
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Figure 5.4 Finite Element Mesh Gradient Near Crack Tip 
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Figure 5.6 Effect of Near Crack Tip Mesh Refinement on 
Crack Tip Response of Homogeneous Beam with 
1.0 Inch Crack 
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Figure 5.11 Dynamic Strain Energy Release Rate 
for Specimen Configuration B 
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Figure 5.16 Dynamic Strain Energy Release Rate 
for Specimen Configuration 6 
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6. SUMMARY 


I 


Dynamic delamination crack propagation in a [90/0]g $ Graphite/ 
Epoxy laminate with an embedded interfacial crack was investigated 
experimentally using high speed photography. The dynamic motion was 
produced by impacting the beam-like laminate specimen with a silicon 
rubber ball. The threshold impact velocities required to initiate 
dynamic crack propagation in laminates with varying initial crack 
positions were determined. The crack propagation speeds were 
estimated from the photographs. 

Experimental results show that the through-the-thickness position 
of the embedded crack can significantly affect the dominant mechanism 
and the threshold impact velocity for the onset of. crack movement. 

If the initial delamination is placed near the top or bottom surface of 
the laminate, local buckling of the delaminated plies may cause 
instability of the crack. If the initial delamination lies on the 
midplane, local buckling does not occur and the initiation of crack 
propagation appears to be dominated by Mode II fracture. The crack 
propagation and arrest observed was seen to be affected by wave 
motion within the delaminated region. 

The contact behavior between the compliant impactor and the 
laminate could not be adequately described by the classical Hertzian 
contact law due to the extent of the deformation and the change in 
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shape undergone by the impactor during the contact interval. An 
empirical model was therefore used to describe the variation of the 
' impact force history with impact velocity. Impact tests of an 
uncracked laminate instrumented with strain gages were used to 
determine the variation of the force history with impact velocity for 
the rubber impactor. It was found that a three-parameter description 
of the force history was necessary to adequately estimate the actual 
contact force. The variation of contact force with time is character- 
ized in this model by two quarter-sine waves of different periods; 
the first represents the loading phase of the force history and the 
second represents the unloading phase. The assumed force history was 
varied in the finite element analysis of the uncracked laminate until 
the calculated strain response sufficiently matched the measured 
response. It was found that the variation of force history with 
impact velocity for the transverse impact of the beam-like laminate 
followed the same trends as in the case of longitudinal impact of a 
prismatic bar. Boundary effects due to wave reflections were not 
significant in the latter case, but were seen to have a considerable 
effect on the force history for the beam specimens. This was due 
primarily to the length of the beam specimens used, and the long 
contact duration of the impactor. As a result, the Force History 
versus Impact Velocity data from these tests did not follow the trends 
predicted by the simple single-degree-of-freedom spring-mass model of 
interaction between the impactor and target, in which wave effects 
are not considered. 
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Calcuation of the dynamic strain energy release rate prior to 
crack extension for specimens of several different geometries indicate 
that there may exist a critical value of this parameter above which 
a delamination crack will become unstable. 

The experimental procedure used here can be improved in several 
ways. Uniformity of the test specimens is necessary whenever the 
results of several tests are to be compared. This was particularly 
important in determining the threshold impact velocities. In fact, 
there were small variations in the length of the composite beam 
specimens, which may have affected the accuracy of the estimate of 
threshold velocity in some cases. More significantly, vibration of 
the long barrel of the air gun caused some variation in the location 
of the impact on the specimens. The magnitude of these variations is 
shown in Tables 3. 1-3.7. The distance between the crack tip and the 
impact site can appreciably affect the impact velocity required to 
cause the crack to extend, so this variation undoubtedly contributed 
to the scatter of the data shown in Figs. 3. 2-3. 7, from which the 
threshold impact velocities were determined. In several cases, the 
accuracy of the estimate of threshold velocity was hindered by the lack 
of sufficient data. The number of different specimen configurations 
tested should be limited, thus allowing more tests of each configuration. 

A more complete evaluation of G c as a criterion for crack 
extension could be performed if the time at which the initial crack 
extension occurred could be determined accurately from the experimental 
data. Comparison of the data and the analysis presented in Chapter 5 
indicate that a higher speed camera would be more suitable for 



experiments of this nature. A camera with a framing rate of one 
frame every ten microseconds would provide data from which a reasonably 
accurate estimate of the time of initial crack extension could be 
obtained. The calculated energy release rate at this time could then 

be used as an estimate of the critical value necessary to extend the 
crack. 

The inverse method used to obtain the parameters describing the 
impact force history may have some inherent inaccuracy. Wave 
reflections from the boundaries of the relatively short impact 
specimens used contributed significantly to the observed strain 
response, and may have obscured that part of the response due solely 
to the contact of the impactor. Therefore, a longer beam would be 
more suitable if this approach is used to obtain impact force 
characteristics of the impactor. 
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